import cv2
import numpy as np
import matplotlib.pyplot as plt

# -------------------------------------巴特沃斯滤波-------------------------------------
def apply_butterworth_filter(image_path, cutoff_frequency, filter_order):
    # 读取图像
    image = cv2.imread(image_path , cv2.IMREAD_GRAYSCALE)
    # 计算图像的中心坐标
    center_x, center_y = image.shape[1] // 2, image.shape[0] // 2

    # 创建巴特沃斯低通滤波器
    butterworth_filter = np.zeros_like(image, dtype=np.float32)
    for y in range(image.shape[0]):
        for x in range(image.shape[1]):
            distance = np.sqrt((x - center_x) ** 2 + (y - center_y) ** 2)
            # 计算巴特沃斯滤波器的值
            butterworth_filter[y, x] = 1 / (1 + (distance / cutoff_frequency) ** (2 * filter_order))
        # cuttoff_frequency是截止频率，filter_order是滤波阶数
        # butterworth_filter 初始化为与输入图像 image 大小相同的全零数组，并设置数据类型为 np.float32。
        # 使用双重循环遍历图像的每个像素 (x, y)。
        # distance 计算每个像素到图像中心的欧氏距离。
        # 每个像素位置的滤波器值按照巴特沃斯滤波器公式计算
        # 其中 D(u, v) 是像素到中心点的距离，D_0 是截止频率，n 是滤波器的阶数。

    # 傅里叶变换图像
    image_frequency = np.fft.fftshift(np.fft.fft2(image))

    # 应用滤波器
    filtered_image_frequency = image_frequency * butterworth_filter
    # 反傅里叶变换
    filtered_image = np.abs(np.fft.ifft2(np.fft.ifftshift(filtered_image_frequency)))

    # 将图像缩放到0-255并转换为uint8类型
    filtered_image = cv2.normalize(filtered_image, None, 0, 255, cv2.NORM_MINMAX)
    filtered_image_uint8 = filtered_image.astype(np.uint8)

    return filtered_image_uint8


def apply_butterworth_highpass_filter(image_path, cutoff_frequency, filter_order):
    # 读取图像
    image = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    # 计算图像的中心坐标
    center_x, center_y = image.shape[1] // 2, image.shape[0] // 2

    # 创建巴特沃斯高通滤波器
    butterworth_filter = np.zeros_like(image, dtype=np.float32)
    for y in range(image.shape[0]):
        for x in range(image.shape[1]):
            distance = np.sqrt((x - center_x) ** 2 + (y - center_y) ** 2)
            # 计算巴特沃斯高通滤波器的值
            butterworth_filter[y, x] = 1 / (1 + (cutoff_frequency / distance) ** (2 * filter_order))
        # cutoff_frequency是截止频率，filter_order是滤波阶数
        # butterworth_filter 初始化为与输入图像 image 大小相同的全零数组，并设置数据类型为 np.float32。
        # 使用双重循环遍历图像的每个像素 (x, y)。
        # distance 计算每个像素到图像中心的欧氏距离。
        # 每个像素位置的滤波器值按照巴特沃斯高通滤波器公式计算
        # 其中 D(u, v) 是像素到中心点的距离，D_0 是截止频率，n 是滤波器的阶数。

    # 傅里叶变换图像
    image_frequency = np.fft.fftshift(np.fft.fft2(image))

    # 应用滤波器
    filtered_image_frequency = image_frequency * butterworth_filter
    # 反傅里叶变换
    filtered_image = np.abs(np.fft.ifft2(np.fft.ifftshift(filtered_image_frequency)))

    # 将图像缩放到0-255并转换为uint8类型
    filtered_image = cv2.normalize(filtered_image, None, 0, 255, cv2.NORM_MINMAX)
    filtered_image_uint8 = filtered_image.astype(np.uint8)

    return filtered_image_uint8


# 主程序
def main():
    # 读取图像
    image_path = 'lena1.png'
    plt.rcParams['font.family'] = 'Times New Roman'

    # 指定截止频率和滤波器阶数
    cutoff_frequency = int(input("输入截止频率："))  # 30
    # 随着截止频率的增加，滤波器的通带宽度也会增加。低截止频率对应于较窄的通带，而高截止频率对应于较宽的通带。
    # 同时，随着截止频率的增加，滤波器的截止斜率也会变平缓。低截止频率下的滤波器更加陡峭，可以更有效地滤除高频成分，但振铃效应也越明显。

    filter_order = int(input("输入阶数："))  # 5
    # 阶数越小，巴特沃斯低通滤波器变化越平滑；阶数越大，其变化越剧烈；若阶数足够高，巴特沃斯低通滤波器应该近似理想低通滤波器，并产生显著的振铃效应。

    # 调用封装函数处理图像
    low_filtered_image = apply_butterworth_filter(image_path, cutoff_frequency, filter_order)
    print(low_filtered_image.shape)
    print(low_filtered_image)
    high_filtered_image = apply_butterworth_highpass_filter(image_path, cutoff_frequency, filter_order)
    # print(high_filtered_image.dtype())
    # 显示滤波后的图像并保存
    plt.figure()
    # 低通滤波
    plt.imshow(low_filtered_image, cmap='gray')
    plt.title(f'Cutoff Frequency: {cutoff_frequency}, Filter Order: {filter_order}')
    plt.xticks([]), plt.yticks([])
    plt.show()
    # # 高通滤波
    # plt.imshow(high_filtered_image, cmap='gray')
    # plt.title(f'Cutoff Frequency: {cutoff_frequency}, Filter Order: {filter_order}')
    # plt.xticks([]), plt.yticks([])
    # plt.show()


# 示例调用
if __name__ == "__main__":
    main()

# -------------------------------------高斯滤波-------------------------------------
import numpy as np
import cv2 as cv


def GLPF(image_path, n, d0):
    """
    Apply Gaussian Low Pass Filter (GLPF) to an image.

    :param image_path: Path to the input image
    :param n: Exponent for shaping the filter
    :param d0: Cutoff frequency
    :return: Filtered image
    """
    # Load and convert the image to grayscale
    image = cv.imread(image_path, cv.IMREAD_GRAYSCALE)
    if image is None:
        raise FileNotFoundError(f"Image at path {image_path} not found.")

    M, N = image.shape

    def generate_glpf(image, d0, n):
        H = np.empty_like(image, dtype=float)
        mid_x = M / 2
        mid_y = N / 2
        for x in range(M):
            for y in range(N):
                d = np.sqrt((x - mid_x) ** 2 + (y - mid_y) ** 2)
                H[x, y] = np.exp(-d ** n / (2 * d0 ** n))
        return H

    def apply_filter(image, filter):
        fft_img = np.fft.fft2(image)  # 空间域转化为频率域，傅里叶二维变换
        fft_img_shift = np.fft.fftshift(fft_img)  # 中心化
        filtered_fft_shift = fft_img_shift * filter
        filtered_fft = np.fft.ifftshift(filtered_fft_shift)  # 反中心化
        filtered_image = np.fft.ifft2(filtered_fft)
        filtered_image = np.abs(filtered_image)
        return filtered_image.astype(np.uint8)

    filter = generate_glpf(image, d0, n)
    filtered_image = apply_filter(image, filter)

    return filtered_image


def GHPF(image_path, n, d0):
    """
    Apply Gaussian High Pass Filter (GHPF) to an image.

    :param image_path: Path to the input image
    :param n: Exponent for shaping the filter
    :param d0: Cutoff frequency
    :return: Filtered image
    """
    # Load and convert the image to grayscale
    image = cv.imread(image_path, cv.IMREAD_GRAYSCALE)
    if image is None:
        raise FileNotFoundError(f"Image at path {image_path} not found.")

    M, N = image.shape

    def generate_ghpf(image, d0, n):
        H = np.empty_like(image, dtype=float)
        mid_x = M / 2
        mid_y = N / 2
        for x in range(M):
            for y in range(N):
                d = np.sqrt((x - mid_x) ** 2 + (y - mid_y) ** 2)
                H[x, y] = 1 - np.exp(-d ** n / (2 * d0 ** n))
        return H

    def apply_filter(image, filter):
        fft_img = np.fft.fft2(image)
        fft_img_shift = np.fft.fftshift(fft_img)
        filtered_fft_shift = fft_img_shift * filter
        filtered_fft = np.fft.ifftshift(filtered_fft_shift)
        filtered_image = np.fft.ifft2(filtered_fft)
        filtered_image = np.abs(filtered_image)
        return filtered_image.astype(np.uint8)

    filter = generate_ghpf(image, d0, n)
    filtered_image = apply_filter(image, filter)

    return filtered_image
# if __name__ == '__main__':
#     image_path = "lena1.png"
#     n = int(input("滤波器的形状的指数:"))
#     # 当n值较大时（例如10），滤波器的衰减速度较快，低频成分受到较强的抑制，因此图像可能会出现更强的模糊效果。
#     d0 = int(input("初始截止频率:"))
#     # 较大的 d0 值会限制低频信息的通过，导致图像变得模糊，因为高频信息被更多地滤除。
#
#     # # 高斯低通滤波
#     # low_filtered_image = GLPF(image_path, n, d0)
#     # print(low_filtered_image.dtype)
#     # cv.imshow('Filtered Image', low_filtered_image)
#     # cv.waitKey(0)
#     # cv.destroyAllWindows()
#
#     # 高斯高通滤波
#     high_filtered_image = GHPF(image_path, n, d0)
#     print(high_filtered_image.dtype)
#     cv.imshow('Filtered Image', high_filtered_image)
#     cv.waitKey(0)
#     cv.destroyAllWindows()
